R Notebook
library(janitor)
##
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(hms)
##
## Attaching package: 'hms'
## The following object is masked from 'package:lubridate':
##
## hms
library(tidyr)
library(stringr)
library(readr)
library(forcats)
library(RcppRoll)
library(tibble)
library(bit64)
## Loading required package: bit
##
## Attaching package: 'bit'
## The following object is masked from 'package:base':
##
## xor
## Attaching package bit64
## package:bit64 (c) 2011-2017 Jens Oehlschlaegel
## creators: integer64 runif64 seq :
## coercion: as.integer64 as.vector as.logical as.integer as.double as.character as.bitstring
## logical operator: ! & | xor != == < <= >= >
## arithmetic operator: + - * / %/% %% ^
## math: sign abs sqrt log log2 log10
## math: floor ceiling trunc round
## querying: is.integer64 is.vector [is.atomic} [length] format print str
## values: is.na is.nan is.finite is.infinite
## aggregation: any all min max range sum prod
## cumulation: diff cummin cummax cumsum cumprod
## access: length<- [ [<- [[ [[<-
## combine: c rep cbind rbind as.data.frame
## WARNING don't use as subscripts
## WARNING semantics differ from integer
## for more help type ?bit64
##
## Attaching package: 'bit64'
## The following objects are masked from 'package:base':
##
## :, %in%, is.double, match, order, rank
library(exploratory)
## Package attached: exploratory v0.3.13. Most recent version available on GitHub: v0.3.16
## You have an OPTION to update the package by typing 'update_exploratory()'. If you do so, make sure to restart R.
##
## Attaching package: 'exploratory'
## The following object is masked from 'package:readr':
##
## read_csv
library(ggplot2)
library(ggridges)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ purrr 0.3.4 ✓ dplyr 1.0.7
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x hms::hms() masks lubridate::hms()
## x lubridate::intersect() masks base::intersect()
## x dplyr::lag() masks stats::lag()
## x exploratory::read_csv() masks readr::read_csv()
## x lubridate::setdiff() masks base::setdiff()
## x lubridate::union() masks base::union()
library(digest)
library(arsenal)
##
## Attaching package: 'arsenal'
## The following object is masked _by_ '.GlobalEnv':
##
## %nin%
## The following object is masked from 'package:lubridate':
##
## is.Date
library(DBI)
library(dplyr)
library(dbplyr)
##
## Attaching package: 'dbplyr'
## The following objects are masked from 'package:dplyr':
##
## ident, sql
library(Hmisc)
## Loading required package: lattice
##
## Attaching package: 'lattice'
## The following object is masked from 'package:exploratory':
##
## histogram
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following object is masked _by_ '.GlobalEnv':
##
## %nin%
## The following object is masked from 'package:arsenal':
##
## %nin%
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
library(knitr)
library(rmarkdown)
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
##
## Attaching package: 'openintro'
## The following object is masked from 'package:survival':
##
## transplant
## The following objects are masked from 'package:lattice':
##
## ethanol, lsegments
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## Registered S3 method overwritten by 'spatstat.geom':
## method from
## print.boxx cli
## spatstat.geom 2.2-2
## Loading required package: spatstat.core
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## spatstat.core 2.3-0
##
## Attaching package: 'spatstat.core'
## The following object is masked from 'package:lattice':
##
## panel.histogram
## The following object is masked from 'package:arsenal':
##
## relrisk
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-0
##
## spatstat 2.2-0 (nickname: 'That's not important right now')
## For an introduction to spatstat, type 'beginner'Purpose
The main purpose of this document is to back the statements made in the manuscript entitled “Variation in Opioid Prescribing Practices Among Female Pelvic Medicine and Reconstructive Surgeons by State from 2013 to 2017: A Retrospective Cohort Study”, but also to describe the methodology applied to the data in order to aggregate and present the results discussed in the manuscript.
Below we set the path to the database to be used as source, as well as each table that will be used throughout the analysis.
db <-'/home/dfcoelho/Downloads/tyler.db'
con <- dbConnect(RSQLite::SQLite(), db)
PUF_tb <- tbl(con, "PUF")
SF_tb <- tbl(con, "SummaryFiles")
FPMRS_tb <- tbl(con, "FPMRS")
ACOG_tb <- tbl(con, "Crosswalk_ACOG_Districts")
DOMS_tb <- tbl(con, "DrugOverdoseMortalityByState")
NPPES_tb <- tbl(con, "NPPES")
OPIOIDS_tb <- tbl(con, "OPIOIDS")Then we get the list of FPMRS practitioners and store in the variable FPMRS, set a list of abbreviations to be filtered out from the list of states. This is followed by a query to build a dataset of opioid prescribing rates at NPI level.
There are a few entries in which bene_count or total_claim_count is equal to NA. That may lead to over count of one of the fields, so we remove any entries in which at least one of the fields is NA.
Here we actually are able to compare the results from the SummaryFiles and PUF tables, which shows that the somehow the datasets don’t actually match after similar aggregation. The only possibility here is to verify the PUF data source.
FPMRS <- FPMRS_tb %>% collect()
## Warning in result_fetch(res@ptr, n = n): Column `Age`: mixed type, first seen
## values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `YearOfBirth`: mixed type, first
## seen values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `LastYearOfResidency`: mixed
## type, first seen values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `ResidencyStartYear`: mixed
## type, first seen values of type integer, coercing other values of type string
FPMRS_NPIs <- unlist(FPMRS$NPI)
excluded_states <- c("AA", "AE", "AP", "AS", "GU", "MP", "VI", "ZZ", "XX")
opioidRX_SF <- SF_tb %>%
filter(NPI %in% FPMRS_NPIs) %>% collect() %>%
filter(State %nin% excluded_states) %>%
mutate_at(vars(OpioidPrescribingRate), parse_number) %>%
mutate_at(vars(Year, NPI, `NPPESProviderLastName`, `NPPESProviderFirstName`,
`Zip`, State, Specialty), factor) %>%
mutate(OpioidClaimCount = as.numeric(OpioidClaimCount)) %>%
mutate(TotalClaimCount = as.numeric(TotalClaimCount)) %>%
mutate(RXRate = OpioidClaimCount/TotalClaimCount) %>%
filter(!is.na(RXRate)) %>% group_by(NPI) %>% mutate(Specialty = "FPMRS") %>%
group_by(NPI) %>%
ungroup() %>% arrange(NPI, Year)
# A list of opioids names (both official and generic)
OPIOIDS_list <- OPIOIDS_tb %>% collect() %>% as.list() %>%
unlist(use.names = FALSE)
opioidRX <-
# get only data from FPMRS using their NPI
PUF_tb %>% filter(NPI %in% FPMRS_NPIs) %>%
# Remove columns that we are not going to use
select(-c(bene_count_ge65, bene_count_ge65_suppress_flag,
total_claim_count_ge65, ge65_suppress_flag,
total_30_day_fill_count_ge65, total_day_supply_ge65,
total_drug_cost_ge65)) %>%
# Remove rows with NAs
filter(!is.na(bene_count)) %>% filter(!is.na(total_claim_count)) %>%
# This forces the data to be collected from the database and be brought to R.
collect() %>%
# Remove entries from with States in the excluded list.
filter(nppes_provider_state %nin% excluded_states) %>%
# Then we rename a few variables
rename(Specialty = specialty_description) %>%
rename(State = nppes_provider_state) %>%
rename(City = nppes_provider_city) %>%
rename(NPPESProviderLastName = nppes_provider_last_org_name) %>%
rename(NPPESProviderFirstName = nppes_provider_first_name) %>%
# We make sure a few variables are factors
mutate_at(vars(Year, NPI, NPPESProviderLastName, NPPESProviderFirstName,
State, Specialty), as.factor) %>%
# Recode values from known factors
mutate(Specialty = recode(
Specialty, `Family Practice` = "Family Medicine",
`Orthopaedic Surgery` = "Orthopedic Surgery",
`Obstetrics/Gynecology` = "Obstetrics & Gynecology")) %>%
mutate(
Specialty = recode(
Specialty, `Obstetrics/Gynecology` = "FPMRS",
Urology = "FPMRS", `Gynecological/Oncology` = "FPMRS",
`Unknown Physician Specialty Code` = "FPMRS",
`Student in an Organized Health Care Education/Training Program` = "FPMRS",
Specialist = "FPMRS", `Osteopathic Manipulative Medicine` = "FPMRS")) %>%
mutate(State = fct_drop(State), Specialty = factor(Specialty)) %>%
# Then we create a variable to hold the opioid count.
# If the drug name or generic name is a match to the drug field, then the
# number held as the total claim count becomes the opioid claim count.
# if the drug is not an opioid, the entries receives a 0 as count
mutate(
opioid_claim_count = ifelse(
drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list,
total_claim_count, 0)) %>%
# Then we group by NPI / Year / Specialty to have an aggregation to a
# physician level
group_by(Year, NPI, NPPESProviderLastName, NPPESProviderFirstName, State,
Specialty) %>%
# And we sum all values for a given physician/year, which show us his stats
# for each period of time
summarise(.groups = "drop",
OpioidClaimCount = sum(opioid_claim_count),
TotalBeneCount = sum(bene_count),
TotalClaimCount = sum(total_claim_count),
RXRate = 100*OpioidClaimCount/TotalClaimCount
) %>% arrange(NPI, Year)
For now, I considered the opioidRX we generated from the PUF files as the correct and continued the analysis.
rmd_total_fpmrs <- opioidRX %>% select(NPI) %>% unique() %>%
ungroup() %>% nrow()
# NPIs of people who prescribed less than ten opioids to medicare beneficiaries
LE10_NPIs <- opioidRX %>% group_by(NPI) %>%
summarise(OpioidClaimCount = round(mean(OpioidClaimCount),0)) %>%
ungroup() %>% filter(OpioidClaimCount < 10)
rmd_fpmrs_data <- opioidRX %>% group_by(NPI) %>%
summarise(OpioidClaimCount = round(mean(OpioidClaimCount),0)) %>%
ungroup() %>% summarise(LE10 = sum(OpioidClaimCount <= 10),
GT10 = sum(OpioidClaimCount > 10)) %>%
mutate(P_LE10 = round(100*LE10/(LE10 + GT10),0),
P_GT10 = round(100*GT10/(LE10 + GT10),0))
fpmrs_GE10 <- opioidRX %>% group_by(NPI) %>%
summarise(OpioidClaimCount = round(mean(OpioidClaimCount), 0)) %>%
ungroup() %>% filter(OpioidClaimCount >= 10) %>% distinct(NPI) %>%
select(NPI) %>% unlist(use.names = FALSE) %>% as.character()A total of 1186 Female Pelvic Medicine and Reconstructive Surgeons (FPMRS) were identified in the Part D Prescriber Public Use File: 646 (54%) prescribed 0-10 opioid claims, and 540 (46%) prescribed more than 10 opioid claims to Medicare Part D patients.
# For the next calculations we need all claims during the period, after dropping
# both unused columns and entries(rows0 with NAs). Dataset include only FPMRS'
claims <- PUF_tb %>%
select(-c(bene_count_ge65, bene_count_ge65_suppress_flag,
total_claim_count_ge65, ge65_suppress_flag,
total_30_day_fill_count_ge65, total_day_supply_ge65,
total_drug_cost_ge65)) %>%
filter(!is.na(bene_count)) %>% filter(!is.na(total_claim_count)) %>%
filter(NPI %in% FPMRS_NPIs) %>% collect()
# We used the claims dataset to calculate the total count of opioids prescribed
total_count <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
summarise(total_bene_count = sum(na.omit(bene_count)),
total_claim_count = sum(na.omit(total_claim_count)))
# then we calculate the basic stats for the claim dataset
stats_claims <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
summarise(
sd_claim = sd(total_claim_count, na.rm = TRUE),
mean_claim = mean(total_claim_count, na.rm = TRUE),
sd_bene = sd(bene_count, na.rm = TRUE),
mean_bene = mean(bene_count, na.rm = TRUE),
mean_supply = mean(total_day_supply /
total_30_day_fill_count, na.rm = TRUE)
)Among FPMRS prescribing at least 10 opioid claims in the Part D Prescriber PUF, a total of 113,201 opioid claims were prescribed to 97,189 beneficiaries (Table 1). Each FPMRS prescribed to a mean (SD) of 32.7 (25.9) opioid claims to a mean of 28.1 (21.1) beneficiaries.
Each beneficiary received a mean of 1.2 opioid claims, with a supply lasting a mean of 5.4 days.
# As we want to present data regarding the statistics of each opioid, we group
# the data by the opioid name and calculate the total claim count and percentage
# of each drug
stats_opioids <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
group_by(drug_name) %>%
summarise(claim_count = sum(total_claim_count), .groups = "drop") %>%
mutate(claim_percentage = round(100 * claim_count / sum(claim_count), 1)) %>%
arrange(desc(claim_count))Combination hydrocodone-acetaminophen accounted for 58,243 opioid claims (52 %), oxycodone acetaminophen 36,457 opioid claims (32%), followed by tramadol hydrochloride 6,750 opioid claims (6 %), and lastly oxycodone 5,116 opioid claims (4%).
# Here the changes concern me a bit, after removing all entries in which
# the beneficiary count or the total claim is NA, the sample size decreased
# significantly and the number of states with mean < 1% near half the number
# of states.
#
# Then I observed what I described on Slack:
# I was checking the states with near zero median and I found out that the values are actually right, since there are a reasonable amount of practitioners that prescribed 0% opioids. In turn, that pushed the median to zero. Closely inspecting the data, I could see that a fair amount of those practitioners had a reduced claimTotal value. So in the median calculation, the practitioner that had 100 claims over the years had the same impact of a practitioner with 10000+ total claim. What if we calculate a weighted median using the total_claim as weight? So the FPMRS that prescribed 10000 drugs (opioids or not) would have a higher impact in the median than a FPMRS what only prescribed 1000 drugs in the same period.
#
# So I calculated the median of OPR but weighted the results by each NPI's
# claim Total, so in a given state the OPR of a professional that prescribed
# hundreds of opioids would impact more than a professional with fewer
# prescriptions. Why do I think that is correct?
# Because medians make more sense to be calculated with a higher granularity level, as if we had one entry per patient and the prescription for that individual. Then, in that situation, if we aggregated to the state level, the median would be "self-weighted". As we don't have this level of information, weighting it by the number of prescriptions would virtually approximate the results.
# opioid prescribing rate = OPR
state_OpioidPR <- claims %>%
mutate(
total_opioid_count = ifelse(
drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list,
total_claim_count, 0)) %>%
group_by(nppes_provider_state, NPI) %>%
summarise(.groups = "drop",
claimTotal = sum(total_claim_count),
opioidTotal = sum(total_opioid_count),
opioidPR = 100*opioidTotal/claimTotal) %>%
group_by(nppes_provider_state) %>%
summarise(
opioidPR = weighted.median(x=opioidPR, w=claimTotal),
claimTotal = sum(claimTotal),
opioidTotal = sum(opioidTotal))
# Then we calculate the top prescribers across all states
top_prescribers <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
group_by(NPI, nppes_provider_state) %>%
summarise(total_claims = sum(total_claim_count), .groups = "drop") %>%
arrange(desc(total_claims))
# And the states with the top prescribers
top_5_states <- top_prescribers %>% head(nrow(.)*0.05) %>%
select(nppes_provider_state) %>% distinct() %>% unlist(use.names = FALSE)
# We select the top three states with higher OPR
top_opr <- state_OpioidPR %>% arrange(desc(opioidPR)) %>%
select(nppes_provider_state) %>% head(3) %>% unlist(use.names = FALSE)
# And the states with OPR lower than 1%. Should we consider get lower values here to narrow down? Perhaps 0.5%? Or just the pure 0%?
bottom_opr <- state_OpioidPR %>% filter(opioidPR < 1) %>%
select(nppes_provider_state) %>% unlist(use.names = FALSE)The median opioid prescribing rate for FPMRS in Alaska, Arkansas, Connecticut, District of Columbia, Delaware, Hawaii, Iowa, Illinois, Kansas, Massachusetts, Maryland, New Hampshire, New York, Pennsylvania, Rhode Island, South Dakota and Vermont is zero, or effectively zero. FPMRS in Mississippi, Maine and Kentucky have the highest median prescribing rates, which are above 6.6%. The top 5% of opioid prescribers were located in Michigan, Mississippi, Texas, Tennessee, West Virginia, Florida, Ohio, Georgia, North Carolina, Missouri, Indiana, Arizona, Nevada, Colorado, California, New Jersey, South Carolina and Alabama.
The median opioid prescription rate for prescriptions covered by Medicare Part D across FPMRS was 1.1% from 2013 to 2017.
A 95% interval for the MEDIAN opioid prescription rate for prescriptions covered by Medicare Part D across FPMRS from 2013 to 2017 was 1.07, 0, 30.72 (Median, Lower, Upper).
A 95% interval for the MEAN opioid prescription rate for prescriptions covered by Medicare Part D across FPMRS from 2013 to 2017 was 5.01, 4.71, 5.3 (Mean, Lower, Upper).
Table:
# calculate the median of claims, escluding claims smaller than 10 opioids
claims_10pp <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(bene_count), " (",
quantile(bene_count) %>% .[c(2, 4)] %>% paste0(collapse = ', '), ")"))
# based on the list of top prescribers, we calculate how many rows is
# approximately 5% of the number of prescribers and we store their NPIs in a
# variable
top_5_prescribers <- top_prescribers %>% head(nrow(.)*0.05) %>% select(NPI) %>%
unlist(use.names = FALSE)
# Then we get the median of prescriptions made by those 5% top prescribers
claims_5pp <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(NPI %in% top_5_prescribers) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(bene_count), " (",
quantile(bene_count) %>% .[c(2, 4)] %>% paste0(collapse = ', '), ")" ))
################################################################################
# Using the statistics for opioids we built before, we create a table with the
# 6 first (most prescribed) and aggregate the rest as others
table_drugs_B_Partial <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(NPI %in% top_5_prescribers) %>%
group_by(drug_name) %>% summarise(claim_count = sum(total_claim_count)) %>%
ungroup() %>%
mutate(claim_percentage = round(100 * claim_count / sum(claim_count), 1)) %>%
arrange(desc(claim_count))
#
table_joined <- stats_opioids %>%
left_join(table_drugs_B_Partial %>%
rename(claim_count_top5=claim_count) %>%
rename(claim_percentage_top5=claim_percentage)
, by="drug_name") %>%
arrange(desc(claim_count))
table_joined_others <- table_joined %>% .[-c(1:6), ] %>%
summarise(
drug_name = "Others", claim_count = sum(claim_count),
claim_count_top5 = sum(claim_count_top5, na.rm = TRUE),
claim_percentage = 0, claim_percentage_top5 = 0)
#
table_1 <- table_joined %>% .[c(1:6), ] %>% rbind(table_joined_others) %>%
summarise(drug_name,
claim_percentage=paste0(claim_count,
" (", round(100*claim_count/sum(claim_count), 2), "%)"),
claim_percentage_top5=paste0(claim_count_top5,
" (", round(100*claim_count_top5/sum(claim_count_top5), 2), "%)")
)
################################################################################
table_total_30_day_fill_count <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_30_day_fill_count,
na.rm = TRUE), "(",
quantile(na.omit(total_30_day_fill_count)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
table_total_30_day_fill_count_5pp <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10,
NPI %in% top_5_prescribers) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_30_day_fill_count,
na.rm = TRUE), "(",
quantile(na.omit(total_30_day_fill_count)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
#-----------------------------------------------------------------------------#
table_total_claim_count <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_claim_count,
na.rm = TRUE), "(",
quantile(na.omit(total_claim_count)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
table_total_claim_count_5pp <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10,
NPI %in% top_5_prescribers) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_claim_count,
na.rm = TRUE), "(",
quantile(na.omit(total_claim_count)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
#-----------------------------------------------------------------------------#
table_total_day_supply <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_day_supply,
na.rm = TRUE), "(",
quantile(na.omit(total_day_supply)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
table_total_day_supply_5pp <- claims %>%
filter(drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list) %>%
filter(total_claim_count > 10,
NPI %in% top_5_prescribers) %>%
summarise(n = n(), `Median (Q1, Q3)` = paste0(median(total_day_supply,
na.rm = TRUE), "(",
quantile(na.omit(total_day_supply)) %>% .[c(2, 4)] %>%
paste0(collapse = ', '), ")"))
The total amount of beneficiaries that received prescriptions from FPMRS (only those with 10+ prescriptions were considered) and from the TOP5% prescribers are:
There are 40 FPMRS that are considered the top 5% prescribers. Also, 785 FPMRS prescribed more than 10 opioids claims.
The top 6 most prescribed drugs are:
total_claim_count > 10 were considered.
total_30_day_fill_count_5pp - Median number of 30 day fill count Opioid Claims prescribed by the top 5% prescribers.
Only prescribers with total_claim_count > 10 were considered.
total_claim_count - Median number of Opioid Claims prescribed
There are only 69 rows in which total_30_day_fill_count != total_claim_count, so total_claim_count_5pp will be virtually equal to table_total_30_day_fill_count_5pp.
total_claim_count_5pp - Median number of Opioid Claims prescribed by the top 5 prescribers.
total_day_supply - Median value of the total amount of days for which a drug was dispensed.
total_day_supply_5pp - Median value of the total amount of days for which a drug was dispensed by the top 5 prescribers.
fpmrs_prescription = PUF_tb %>%
filter(!is.na(bene_count)) %>% filter(!is.na(total_claim_count)) %>%
filter(NPI %in% FPMRS_NPIs) %>% collect() %>%
mutate(
total_opioid_count = ifelse(
drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list,
total_claim_count, 0)) %>%
summarise(
total = sum(total_claim_count, na.rm = TRUE),
total_opioid = sum(total_opioid_count, na.rm = TRUE)
) %>% mutate(prescription_rate = 100*(total_opioid/total))
The opioid prescription rate across FPMRS, whose prescriptions were covered by Medicare Part D, from 2013 to 2017 was 3.6%.
total_opioid_prescription = PUF_tb %>%
filter(!is.na(bene_count)) %>% filter(!is.na(total_claim_count)) %>%
mutate(
total_opioid_count = ifelse(
drug_name %in% OPIOIDS_list | generic_name %in% OPIOIDS_list,
total_claim_count, 0)) %>%
summarise(
total = sum(total_claim_count, na.rm = TRUE),
total_opioid = sum(total_opioid_count, na.rm = TRUE)
) %>% mutate(prescription_rate = 100*(total_opioid/total)) %>% collect()These FPMRS providers prescribed a total of 113,201 opioid prescriptions or 0.038% of all opioid prescriptions covered by Medicare Part D from 2013 to 2017.
source('Ridge_and_HeatMap_DB_Rev1.R')
## Warning in result_fetch(res@ptr, n = n): Column `Age`: mixed type, first seen
## values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `YearOfBirth`: mixed type, first
## seen values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `LastYearOfResidency`: mixed
## type, first seen values of type integer, coercing other values of type string
## Warning in result_fetch(res@ptr, n = n): Column `ResidencyStartYear`: mixed
## type, first seen values of type integer, coercing other values of type string
## `summarise()` has grouped output by 'NPI', 'State'. You can override using the `.groups` argument.
## Picking joint bandwidth of 0.0486
## `summarise()` has grouped output by 'State'. You can override using the `.groups` argument.
# The command below will convert TIFF to PNG (requires imagemagick)
system("sh Tiff2Png.sh")Ridges were not included for states with less than 4 FPMRS listed in the PUF file.
